1. LECTURA DE DADES

Un cop descarregades les dades, les llegim i les grafiquem en forma de sèrie temporal. Visualitzam la sèrie de Mallorca, Menorca, Eivissa i Formentera, i les Illes Balears en una mateixa gràfica:

turisme <- read_excel("IBESTAT.xls")
turisme$Data <- gsub("M","-",turisme$Data)
gastos.ts<-ts(turisme[-1], start=c(2015,10), frequency = 12)
plot(gastos.ts, plot.type = "single", col = 1:ncol(gastos.ts), ylim=c(0,1300), xlab="Any", ylab="Despeses mensuals en €", main='Series de cada illa')
legend("bottomright", colnames(gastos.ts), col=1:ncol(gastos.ts), lty=1.5, cex=.85)

Les dades van des d’octubre de 2015 a novembre de 2022 (llavors tenim 6 cicles complets). Podem observar que totes les illes segueixen més o menys un mateix patró (amb una diferència entre les series de Menorca i Eivissa i Formentera davant les de Mallorca i les Illes Balears) i l’efecte de la COVID és veu clarament en 2020. També observam una clara estacionalitat que pareix tenir els seus valors més alts en estiu i els més baixos en hivern.

En aquest document estudiarem l’illa de Mallorca, diguem mall al conjunt de dades que fan referència a l’illa de Mallorca.

mall <- data.frame(x=1:86, y=turisme$Mall) 
head(mall)
##   x      y
## 1 1 895.45
## 2 2 850.73
## 3 3 834.39
## 4 4 773.57
## 5 5 858.62
## 6 6 814.16

La variable \(X\) ens indica el mes en que ens trobam (hem de tenir en compte que les dades comencen el 10 de novembre 2015, llavors el valor 1 correspon a aquest més i segueix fins al mes 86, que correspon a novembre de 2022), mentre que la variable \(Y\) ens indica les despeses dels turistes en aquells mes en euros.

Abans d’aplicar algun model i posar-nos a analitzar les dades, visualitzem-les per tenir una idea de qué ens podem trobar a l’hora d’aplicar els models:

serie_mall <- ts(mall$y,start=c(2015,10),frequency = 12)
plot.ts(serie_mall, main="Illa de Mallorca", ylab="Despeses mensuals en €", xlab="Any")

Per veure millor l’estacionalitat de la sèrie visualitzem cada període mensualment:

seasonplot(serie_mall, col=c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),year.labels=TRUE, main="Estacionalitat de Mallorca", xlab="Mes", ylab=" Despeses en €")

legend("bottomright",
       legend = c(2015, 2016, 2017, 2018, 2019, 2020, 2021, 2022),
       fill =  c("brown", "blue","red", "orange", "pink", "purple", "yellow","green"),      
       border = "black") 

Observam que efectivament, així com es comenta en el treball, el període de confinament va ser de principi de març, on comença a decréixer el valor dels preus, fins a finals de juny, on pareix que s’han tornat a recuperar els valors anteriors a la COVID. A més, també observam que les temporades d’estiu són més altes que les de l’hivern.

Dividirem la sèrie en tres trams:

  • Des de 10-2015 fins 12-2018 (serie1_mall), amb el que predirem un període. És a dir, predirem de 1-2019 a 12-2019. (amb aquest comprovam que els mètodes funcionen)

  • Des de 10-2015 fins 12-2019 (serie2_mall), amb la que predirem el període de la COVID (és a dir, l’any 2020)

  • Des de 10-2015 fins 12-2020 (serie3_mall), amb la que predirem el període després de la COVID (és a dir, l’any 2021)

Visualitzem-los:

serie1_mall <- ts(mall$y[1:39],start=c(2015,10),frequency = 12)
serie2_mall <- ts(mall$y[1:51],start=c(2015,10),frequency = 12)
serie3_mall <- ts(mall$y[1:63],start=c(2015,10),frequency = 12)

par(mfrow=c(2,2), mar=c(4,4,4,1)+.1)
plot.ts(serie1_mall, main="Sèrie abans de la COVID \nmenys un cicle" , xlab="Any", ylab="despeses mensuals")
plot.ts(serie2_mall, main="Sèrie abans de la COVID", xlab="Any", ylab="despeses mensuals")
plot.ts(serie3_mall, main="Sèrie abans i durant de la COVID", xlab="Any", ylab="despeses mensuals")
plot.ts(serie_mall, main="Sèrie completa", xlab="Any", ylab="despeses mensuals")


2. AJUST DELS MODELS:

Abans d’aplicar algun model, estudiem un poc el tram a analitzar:

serie_mall.lm <- lm(y~x, data=mall) 

plot.ts(mall$y, main = "Mallorca ", ylab = "Despeses menusals en €", xlab="Índex de cada mes")

#dibuixam la recta de regressió sobre els punts
abline(serie_mall.lm, col='red')

Si dibuixam la recta de regressió sobre les nostres dades, tot i que aquestes estan molt disperses i no s’ajusten bé a la recta, podem observar que la tendència és una mica creixent, encara que es manté més o menys constant.

boxplot(serie_mall~cycle(serie_mall), xlab = "Mes", ylab = "Despeses en €", main="Boxplot de Mallorca")

Podem observar també la presència d’estacionalitat, que prenen els seus valors màxims a la temporada d’estiu i els seus mínims en l’hivern, fet que correspon amb les temporades turístiques a l’illa.

Aplicarem diversos models per ajustar la nostra sèrie i fer-ne una predicció per llavors comparar quin és el millor.

Vegem ara com actua cada model:


2.1. ABANS DE LA COVID


- REGRESSIÓ SEGMENTADA

Com hem comentat anteriorment la recta de regressió no s’ajusta bé a les dades, de fet el valor del \(R^2\) és molt baix:

serie1_mall.lm <- lm(y~x, data=mall[1:39,])
summary(serie1_mall.lm)$adj.r.squared
## [1] 0.02188284

Per això utilitzam el paquet segmented per ajustar a una recta de regressió a trossos.

Anem a utilitzar la funció selgmented() per veure quants de punts de canvi detecta:

punts_canvi_serie1_mall <-selgmented(serie1_mall.lm, Kmax = 10, type="bic")
## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 .. 
## 
## BIC to detect no. of breakpoints:
##        0        1        2        3        4        5        6        7 
## 370.5453 374.5279 381.6098 376.5556 372.6382 357.6041 351.5419 337.9820 
##        8        9       10 
## 344.0487 349.0722 355.4233 
## 
## No. of selected breakpoints: 7

Aplicam la funció segmented() que ens calcula la regressió segmentada:

serie1_mall.seg <- segmented(serie1_mall.lm, seg.Z=~x, psi = c(3,10,17,23,28,35))

Aquesta funció ens demana que introduïm els valors on es troben els punts de canvi, i aquesta ens dóna el valor estimat. Aquests punts de canvi són:

summary(serie1_mall.seg)$psi
##        Initial      Est.    St.Err
## psi1.x       3  6.557347 0.7808007
## psi2.x      10 10.579679 0.4882573
## psi3.x      17 16.351485 0.4500235
## psi4.x      23 22.934881 0.4249014
## psi5.x      28 27.781744 0.4721648
## psi6.x      35 34.980655 0.4401346

Que corresponen, aproximadament, a: 4-2016, 8-2016, 1-2017, 8-2017, 1-2018, 8-2018. Llavors tenim que els punts més alts es donen per agost i els més baixos per gener.

Obtenim que el valor de \(R^2\) per a la regressió segmentada és bastant alt

summary(serie1_mall.seg)$adj.r.squared
## [1] 0.8063167

Anem a visualitzar la regressió segmentada damunt les nostres dades

#Per graficar-ho

p_serie1_mall <- ggplot(mall[1:39,], aes(x=x, y=y), color="black") + 
  geom_line()+ #dibuixam les nostres dades en línia contínua
  ggtitle('Regressió lineal i segmentada sobre les dades de Mallorca \nabans de la COVID') + 
  xlab('Índex del mes') + 
  ylab('Despeses mensuals en €') +
  ylim(c(0,1300))
  

my.coef1_mall <- coef(serie1_mall.lm) #coeficients de la recta de regressió lineal

p_serie1_mall <- p_serie1_mall + geom_abline(intercept=my.coef1_mall[1], slope=my.coef1_mall[2], color="green") #afegim la recta de regressió lineal

my.model1_mall <- data.frame(x=1:39, y=fitted(serie1_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada

p_serie1_mall <- p_serie1_mall + geom_line(data=my.model1_mall, aes(x=x,y=y), color="red") #afegim la recta de regressió segmentada

my.lines1_mall <- serie1_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats

p_serie1_mall <- p_serie1_mall + geom_vline(xintercept=my.lines1_mall, linetype="dashed") #línies verticals en els punts de canvi

p_serie1_mall <- p_serie1_mall + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

ggplotly(p_serie1_mall)

Visualment també es pot observar que la recta de regressió segmentada s’ajusta millor a les nostres dades.

Anem a calcular les equacions d’aquestes rectes, sabem que les rectes tenen la forma \(y=mx+n\), on \(m\) és la pendent i \(n\) el valor de tall de l’eix y.

Hi ha una funció del paquet segmented que ens dóna aquests valors de \(n\):

intercept(serie1_mall.seg)
## $x
##               Est.
## intercept1  882.18
## intercept2  354.63
## intercept3 1740.80
## intercept4 -300.59
## intercept5 2844.30
## intercept6 -606.28
## intercept7 3780.00

També tenim una altra funció que ens calcula les pendents:

slope(serie1_mall.seg)
## $x
##           Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -12.674 11.4980 -1.1023   -36.355    11.006
## slope2  67.777 21.5110  3.1509    23.475   112.080
## slope3 -63.248 11.4980 -5.5008   -86.928   -39.567
## slope4  61.599 11.4980  5.3574    37.919    85.279
## slope5 -75.525 15.2100 -4.9654  -106.850   -44.199
## slope6  48.679  9.0899  5.3553    29.958    67.400
## slope7 -76.713 15.2100 -5.0435  -108.040   -45.387

Aleshores, la nostra recta segmentada és:

\(\hat{y}= \left\{ \begin{array}{lcc} -12.674x + 882.18 & si & x \leq \psi_1 \\ \\ 67.777x + 354.63 & si & \psi_1 < x \leq \psi_2 \\ \\ -63.248x + 1740.8 & si & \psi_2 < x \leq \psi_3 \\ \\ 61.599x - 300.59 & si & \psi_3 < x \leq \psi_4 \\ \\ -75.525x + 2844.3 & si & \psi_4 < x \leq \psi_5 \\ \\ 48.679x - 606.28 & si & \psi_5 < x \leq \psi_6 \\ \\ -76.713x + 3780 & si & \psi_6 < x \\ \end{array} \right.\)

on \(\psi_1= 6.56\), \(\psi_2= 10.58\), \(\psi_3= 16.35\), \(\psi_4= 22.93\), \(\psi_5= 27.78\) i \(\psi_6= 34.98\)

Vegem els errors del model (ens interessa el RMSE):

accuracy(serie1_mall.seg)
##                        ME     RMSE     MAE        MPE     MAPE      MASE
## Training set 2.915581e-15 38.51009 32.1737 -0.1792044 3.633043 0.3559833

Anem a fer la predicció de 1-2019 a 12-2019.

El darrer punt de canvi que tenim és en agost de 2018 i sabem, seguint el patró obtingut amb les dades d’entrenament, que els següents és donaran en gener de 2019, agost de 2019 i gener de 2020. Necessitam calcular les pendents de les rectes d’entre agost de 2018 i gener de 2020, per poder fer la predicció i calcular l’error.

Per calcular les pendents de les noves rectes farem la mitjana de les pendents anteriors. De les pendents ja calculades en el model obviarem la primera i la darrera, ja que no són vàlides. Per tant,

  • La pendent de 8-2018 a 1-2019 i de 8-2019 a 1-2020 serà : (-63.248 - 75.525)/2 = -69.387
  • La pendent de 1-2019 a 8-2019 serà : (66.777 + 61.599 + 48.679)/3 =59.018

Llavors els nous segments predits seran

\(\hat{y}= \left\{ \begin{array}{lcc} -69.387x + n_1 & si & x \leq \psi_7 \\ \\ 59.018x + n_2 & si & \psi_7 < x \leq \psi_8 \\ \\ -69.387x + n_3 & si & \psi_8 < x \\ \end{array} \right.\)

on \(\psi_7 = 40\) i \(\psi_8 = 47\).

Calculem aquests \(n_i\), \(i\in \{1,2,3\}\):

La darrera recta vàlida de la regressió segmentada és \(y = 48.679x - 606.28\)

Per \(n_1\);

Si \(x = 34.98\), que és el darrer punt de canvi, la \(y = 48.679 \cdot 34.98 - 606.28 = 1096.511\).Per tant, el punt d’intersecció de la primera recta a predir és \(1096.511 = -69.387 \cdot 34.98 + n_1\). Llavors \(n_1 = 3523.67\).

Per \(n_2\);

Si \(x= 40\), que és el primer punt de canvi predit, tenim \(y = -69.387 \cdot 40 + 3523.67 = 748.19\). Per tant, el segon punt d’intersecció amb l’eix OY és \(748.19 = 59.018 \cdot 40 + n_2\), \(n_2 = -1612.53\).

Per \(n_3\);

Si \(x=47\), que és el segon punt de canvi predit, \(y = 59.018 \cdot 47 -1612.53 = 1161.315\). Per tant, el tercer punt d’intersecció amb l’eix OY és \(1161.315 = -69.387 \cdot 47 + n_3\), \(n_3 = 4422.48\).

Aleshores, la predicció per al 2019 amb la regressió segmentada és

\(\hat{y}= \left\{ \begin{array}{lcc} -69.387x + 3523.67 & si & x \leq \psi_7 \\ \\ 59.018x -1612.53 & si & \psi_7 < x \leq \psi_8 \\ \\ -69.387x + 4422.48 & si & \psi_8 < x \\ \end{array} \right.\)

on \(\psi_7 = 40\) i \(\psi_8 = 47\).

#Per graficar-ho

p1_serie_mall <- ggplot(mall[1:51,], aes(x=x, y=y)) + geom_line()+
  ylim(c(500,1300)) +#dibuixam les nostres dades en línia contínua
  ggtitle('Predicció de Mallorca amb el model \nde regressió segmentada abans de la COVID') + 
  xlab('Índex del mes') + 
  ylab('Despeses mensuals en €')

my.model1_mall <- data.frame(x=1:39, y=fitted(serie1_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada

p1_serie_mall <- p1_serie_mall + geom_line(data=my.model1_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines1_mall <- serie1_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats

p1_serie_mall <- p1_serie_mall + geom_vline(xintercept=my.lines1_mall, linetype="dashed") #línies verticals en els punts de canvi

p1_serie_mall <- p1_serie_mall+ geom_vline(xintercept=0) + geom_hline(yintercept=500) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

p1_serie_mall <- p1_serie_mall + geom_abline(intercept = 3523.67, slope=-69.387, colour="green") +
                geom_abline(intercept = -1612.53, slope=59.018, colour="blue") + 
                geom_abline(intercept = 4422.48, slope=-69.387, colour="orange")
  

ggplotly(p1_serie_mall)

Observam que les noves rectes predites s’ajusten més o menys bé a la sèrie.

Calculem l’error de la predicció:

o1_mall<-c(serie_mall[40:51]) #observacions reals de la sèrie  de gener de 2019 a desembre de 2019

v1=c(40:47)
f1 <- sapply(v1, function(x) 59.018*x-1612.53) #predicció de gener 2019 a agost 2019  

v2=c(48:51)
f2 <- sapply(v2, function(x) -69.387*x + 4422.48)#predicció de setembre 2019 a desembre 2019 

p1_mall<- c(f1,f2) #predicció de gener de 2019 a desembre de 2019

sqrt(sum((p1_mall-o1_mall)^2)/12) #error de la predicció
## [1] 52.35583

Aleshores aquest model pareix ser bastant bo, ja que per fer la predicció un any endavant només tenim un error d’uns 43 euros.


- DECOMPOSE

Recordem la sèrie

plot.ts(serie1_mall, main= "Mallorca abans de la COVID", xlab="Any", ylab="Despeses mensuals en €")

Ja hem dit anteriorment que té una tendència mínima i podem observar, també, que no hi ha variabilitat. Llavors podem suposar un model additiu.

El mètode clàssic de descomposició separa la sèrie en subseries corresponents a la tendència, la estacionalitat i el renou.

Aquestes components es poden combinar de manera additiva o multiplicativa. En el nostre cas utilitzam el model additiu: \(y_t = \mu_t + S_t + a_t\)

decompose_s1_mall<-decompose(serie1_mall)
plot(decompose_s1_mall, xlab="Any")

El decompose, per calcular aquestes noves tendències, el que fa és agafar les sis tendències anteriors i les sis posteriors de la sèrie original i en fa la mitjana. És per això que la primera que obtenim és en abril de 2016 i la darrera en juny de 2018. Notem que per a la predicció ens quedarem amb el valor de la darrera tendència del decompose.

t1_mall <- decompose_s1_mall$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t1_mall
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2015                                                                        
## 2016       NA       NA       NA 899.8458 899.0708 896.9637 895.5496 888.9975
## 2017 891.8125 896.2783 901.1758 905.1896 907.7217 908.4975 909.7354 914.6108
## 2018 924.3500 924.6162 923.1429 921.8554 921.3479 920.0171       NA       NA
##           Sep      Oct      Nov      Dec
## 2015                NA       NA       NA
## 2016 885.1633 886.3112 887.5750 888.9183
## 2017 920.9775 923.1212 922.9363 923.5758
## 2018       NA       NA       NA       NA

Els valors de les components estacionals els calcula fent la mitjana per mesos, és a dir, per calcular la componen de gener, agafa els valors de tots els geners que tenim a la sèrie original i en fa la mitjana. Per tant, només tenim 12 valors, un per a cada mes.

s1_mall <- decompose_s1_mall$seasonal #components estacionals
s1_mall <-  s1_mall[4:15] #estacionalitat de gener a desembre
s1_mall
##  [1] -174.312598 -107.883640  -54.675723  -66.333293    6.116846   47.367541
##  [7]  169.981152  195.059485   82.053235   67.382402 -116.501973  -48.253432

Anem a fer la predicció d’aquest model

a1_mall <- c(s1_mall[7:12],s1_mall) #estacionalitat de juliol 2018 a desembre 2019 (és per poder fer la predicció)

pred1_decompose_mall <- sapply(a1_mall, function(x) 920.0171 + x) #predicció de juliol  de 2018 a desembre 2019 (el valor 920.0171 és la darrera tendència obtinguda amb el decompose)
p1_dec_mall<-c(serie1_mall[1:33], pred1_decompose_mall) #valors de la sèrie original més els predits

A1_mall<- data.frame("x" = mall[1:51,]$x, "y" = mall[1:51,]$y, "p"= p1_dec_mall) #cream un data frame on hi trobam una columna amb els mesos, x, una amb els valors originals de la sèrie, y, i una altra amb les prediccions, p.

#ho dibuxam
grafica1_mall_dec <- ggplot(A1_mall)+ 
  ylim(c(600,1200)) + 
  geom_line(aes(x,p), color="red")+
  geom_line(aes(x,y))+ geom_vline(xintercept=0) +
  geom_hline(yintercept=600)+ 
  labs(title="Predicció de Mallorca abans de la COVID amb el model de descomposició", x="Índex del mes", y="Despeses mensuals en €")+ 
  theme(panel.background = element_blank())

grafica1_mall_dec

Podem observar que aquest model fa una predicció bastant bona. No obstant això observam que, com el quart cicle és un poc més alt que els anteriors, la predicció no arriba a captar els valors més alts del darrer cicle.

Calculem l’error de la predicció:

(Notem que la predicció és des de juliol de 2018 a desembre de 2019 i per calcular l’error només volem de gener de 2019 a desembre de 2019.)

sqrt(sum((c(serie_mall[40:51]- pred1_decompose_mall[7:18]))^2)/12) #error predicció de gener a desembre de 2019
## [1] 32.40833

Aleshores l’error que és comet amb el model clàssic de descomposició és de 32 euros, aproximadament.


Notem que també tenim una altra instrucció a R per fer prediccions d’una sèrie, la funció predict(). Aquesta és basa en un model ETS, anem a veure la predicció:

prediccio1_mall <- predict(serie1_mall,h=12)
plot(prediccio1_mall, xlab="Any", ylab="Despeses mensuals en €")

Pareix que la predicció és bastant bona, ja que el cicle predit segueix un mateix patró que els anteriors. Anem a veure la predicció juntament amb la sèrie original:

df_prediccio1_mall <- data.frame(prediccio1_mall)
p1_ets_mall <- data.frame("x"= 1:51, "PointForecast"=c(serie1_mall,df_prediccio1_mall$Point.Forecast), "Lo80" = c(rep(NA,39),df_prediccio1_mall$Lo.80), "Hi80" = c(rep(NA,39),df_prediccio1_mall$Hi.80), "Lo95" = c(rep(NA,39),df_prediccio1_mall$Lo.95),"Hi95" = c(rep(NA,39),df_prediccio1_mall$Hi.95))

grafica_pred1_ets <- ggplot((mall[1:51,]))+ ylim(c(600,1200))+
  geom_ribbon(data = p1_ets_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data =p1_ets_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = p1_ets_mall, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció del model ETS(A,N,A) a Mallorca abans de ls COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=600)+ theme(panel.background = element_blank())

grafica_pred1_ets

Podem observar que la predicció és bastant bona, ja que continua seguint un mateix patró. I l’error comés seria de 37 euros.

accuracy(prediccio1_mall,serie2_mall)
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set  1.886318 26.85308 22.27401 0.1691627 2.499489 0.6617645
## Test set     19.065194 37.02366 31.47709 1.6792842 3.289203 0.9351895
##                     ACF1 Theil's U
## Training set -0.01565369        NA
## Test set      0.43720293 0.4375391


- SARIMA

Anem a veure quin model obtenim considerant un model estacional.

Pel que hem vist anteriorment, podem considerar que no hi ha tendència, llavors no fa falta fer cap diferència a la part regular, no obstant, sí que cal fer una diferenciació d’orde 12 per l’estacionalitat. Vegem l’ACF i el PACF de la sèrie a predir:

par(mfrow=c(1,2))
acf(serie1_mall)
pacf(serie1_mall)

L’ACF ens està indicant un model de mesures mòbils d’ordre 2, MA(2). És a dir, té en compte els errors comesos en els dos mesos anteriors a l’hora de fer la predicció. A més, el PACF ens indica un model autoregressiu d’ordre 1, AR(1). Aquest fet ens està indicant que a l’hora de fer la predicció es té en compte el valor obtingut al mes anterior.

Llavors, per a la part regular obtenim un ARIMA(1,0,2).

Fem una diferenciació a la part estacional, és a dir, d’ordre 12

serie1_mall_dif <- diff(serie1_mall,12)
plot(serie1_mall_dif, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")

Aquesta és la sèrie sense estacionalitat, vegem com es modifiquen l’ACF i el PACF.

par(mfrow=c(2,1))
acf(serie1_mall_dif, lag=36)
pacf(serie1_mall_dif,lag=36)

Observam que no hi ha barres significatives, llavors no tenim cap model autoregressiu ni de mesures mòbils a la part estacional. És a dir, no es tenen en compte els valors obtinguts dels cicles anteriors.

Per tant hem obtingut un ARIMA(1,0,2)(0,1,0)[12].

Vegem quin model els proposa R:

auto.arima(serie1_mall)
## Series: serie1_mall 
## ARIMA(0,0,0)(1,1,0)[12] with drift 
## 
## Coefficients:
##          sar1   drift
##       -0.5576  0.8407
## s.e.   0.1924  0.4397
## 
## sigma^2 = 1342:  log likelihood = -136.72
## AIC=279.45   AICc=280.49   BIC=283.33

R ens proposa un ARIMA(0,0,0)(1,1,0)[12], és a dir només fa una diferenciació a la part estacional i considera que per predir un mes s’han de tenir en compte els errors d’aquest mateix mes però en el cicle anterior.

Per tant les propostes de model ARIMA són:

model1_mall<-arima(serie1_mall, order=c(1,0,2), seasonal = list(order=c(0,1,0), period=12))

model2_mall <- arima(serie1_mall, order=c(0,0,0), seasonal = list( order=c(1,1,0), period=12))

Els errors d’aquests dos models obtinguts són

accuracy(model1_mall)
##                    ME     RMSE      MAE      MPE     MAPE      MASE        ACF1
## Training set 4.916634 32.66904 22.17796 0.380845 2.549795 0.3062364 -0.03442627
accuracy(model2_mall)
##                    ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set 8.180359 32.10763 21.18521 0.6723209 2.372129 0.2925283 0.2063501

Observam que no hi ha molta diferència d’un model a un altre.

Anem a fer la predicció de 2019 d’aquests models :

#ens torna un data frame amb els punts predits, un interval del 80% i un altre del 95%
fc_m1_mall <-forecast(model1_mall, h=12) 
fc_m2_mall <- forecast(model2_mall,h=12)

Visualitzem-les:

#MODEL 1

pre1 <- data.frame("Point Forecast" = serie1_mall, "Lo 80" = rep(NA,39), "Hi 80"=  rep(NA,39), "Lo 95" =  rep(NA,39), "Hi 95" =  rep(NA,39)) #dades abans de la predicció amb NA als intervals ja que només els volem per la predicció
 
pred_m1_mall <-data.frame(fc_m1_mall) 

#vectors amb les dades de la sèrie i la predicció amb els intervals
a1_mall <- c(pre1$Point.Forecast,pred_m1_mall$Point.Forecast) 
b1_mall <- c(pre1$Lo.80, pred_m1_mall$Lo.80)
c1_mall <- c(pre1$Hi.80, pred_m1_mall$Hi.80)
d1_mall <- c(pre1$Lo.95, pred_m1_mall$Lo.95)
e1_mall <- c(pre1$Hi.95, pred_m1_mall$Hi.95)

grafica_m1_mall <- data.frame("x" = 1:51, "PointForecast" = a1_mall, "Lo80" = b1_mall, "Hi80"= c1_mall, "Lo95" = d1_mall, "Hi95" = e1_mall)

grafica1_mall <- ggplot(mall[1:51,])+ ylim(c(600,1200))+
  geom_ribbon(data = grafica_m1_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m1_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m1_mall, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció de Mallorca amb el model ARIMA (1,0,2)(0,1,0)[12] abans de la COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=600)+ theme(panel.background = element_blank())
grafica1_mall

Podem observar que la predicció (línia blava) és bastant bona ja que és molt propera a la de les dades reals (la vermella). A més, les bandes dels intervals de confiança no són molt amplis, el que ens indica que la probabilitat d’error no és molt gran.

#MODEL 2

pred_m2_mall <-data.frame(fc_m2_mall)

a2_mall <- c(pre1$Point.Forecast,pred_m2_mall$Point.Forecast)
b2_mall <- c(pre1$Lo.80, pred_m2_mall$Lo.80)
c2_mall <- c(pre1$Hi.80, pred_m2_mall$Hi.80)
d2_mall <- c(pre1$Lo.95, pred_m2_mall$Lo.95)
e2_mall <- c(pre1$Hi.95, pred_m2_mall$Hi.95)

grafica_m2_mall <- data.frame("x" = 1:51, "PointForecast" = a2_mall, "Lo80" = b2_mall, "Hi80"= c2_mall, "Lo95" = d2_mall, "Hi95" = e2_mall)


grafica2_mall <- ggplot(mall[1:51,])+ ylim(c(600,1200))+
  geom_ribbon(data = grafica_m2_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m2_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m2_mall, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció de Mallorca amb el model ARIMA (0,0,0)(1,1,0)[12] abans de la COVID", x="índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=600)+ theme(panel.background = element_blank())

grafica2_mall

Amb aquest model la predicció pareix un poc més bona, ja que les bandes blaves són més fines i la línia blava, que són els valors predits, s’ajusten millor als reals (els vermells).

Vegem quins són els errors de la predicció:

accuracy(fc_m1_mall, serie_mall[40:51], h=12)
##                     ME     RMSE      MAE      MPE     MAPE      MASE
## Training set  4.916634 32.66904 22.17796 0.380845 2.549795 0.3062364
## Test set     13.314448 39.50200 32.06322 1.095781 3.483925 0.4427334
##                     ACF1
## Training set -0.03442627
## Test set              NA
accuracy(fc_m2_mall,serie_mall[40:51], h=12)
##                     ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set  8.180359 32.10763 21.18521 0.6723209 2.372129 0.2925283 0.2063501
## Test set     17.204294 30.21333 24.52273 1.6109194 2.569169 0.3386133        NA

Efectivament, com hem vist que pareixia amb les gràfiques, la predicció del segon model dona menys error (30.21) que el primer(39.5), però la diferència és mínima.

Els errors o residus són útils per comprovar si un model ha capturat adequadament la informació de les dades. Un bon mètode de pronòstic produirà residus amb les següents propietats:

  1. Els residus no estan correlacionats. Si hi ha correlació, queda estructura per eliminar.

  2. Els residus tenen mitjana zero.

  3. La variància és constant.

  4. Es distribueixen normalment.

Estudiem-los:

checkresiduals(fc_m1_mall)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,2)(0,1,0)[12]
## Q* = 3.4598, df = 5, p-value = 0.6295
## 
## Model df: 3.   Total lags used: 8
checkresiduals(fc_m2_mall)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0)(1,1,0)[12]
## Q* = 6.5782, df = 7, p-value = 0.4741
## 
## Model df: 1.   Total lags used: 8

Amb les gràfiques del checkresiduals podem observar que, en els dos models, la sèrie dels errors no mostra cap estructura i les barres dels ACF es troben dins les bandes de confiança. A més, s’observa que les distribucions dels errors no pareixen seguir una normal. Aquesta funció també ens proporciona el p-valor del test de Ljung-Box, que avalua si els residus del model s’ajusten a un procés de renou blanc, és a dir, si els residus són independents i no correlacionats. El model ARIMA(1,0,2)(0,1,0)[12] pareix seguir millor l’estructura de la gaussiana.

Vegem els qqPlot i apliquem tests de normalitat:

par(mfrow=c(1,2))
qqPlot(fc_m1_mall$residuals, id=FALSE, mean=mean(fc_m1_mall$residuals), sd=sd(fc_m1_mall$residuals))
qqPlot(fc_m2_mall$residuals, id=FALSE, mean=mean(fc_m2_mall$residuals), sd=sd(fc_m2_mall$residuals))

Podem observar que en el primer model tenim més valors centrats al 0 i menys valors fora de la banda. Tot i així els valors no es troben sobre la diagonal, fet que ens indica que no hi ha normalitat.

shapiro.test(fc_m1_mall$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m1_mall$residuals
## W = 0.91601, p-value = 0.006562
shapiro.test(fc_m2_mall$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m2_mall$residuals
## W = 0.89116, p-value = 0.001243
lillie.test(fc_m1_mall$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m1_mall$residuals
## D = 0.17804, p-value = 0.003142
lillie.test(fc_m2_mall$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m2_mall$residuals
## D = 0.25238, p-value = 1.126e-06

Tots els tests ens donen un valor molt baix, concretament menor a 0,05. Fet que ens indica que no hem d’acceptar la hipòtesi de que la distribució segueix una normal.


Vegem un resum dels errors que cometen cada un dels models anteriors:

reg. segmentada descomposició clàssica ETS(A,N,A) ARIMA (1,0,2)(0,1,0)[12] ARIMA (0,0,0)(1,1,0)[12]
error model 38.51 26.85 32.6690 32.1076
error predicció 52.356 32.40833 37.02366 39.5020 30.21333


2.2. DURANT LA COVID


- REGRESSIÓ SEGMENTADA

Ara el valor de \(R^2\) per a la regressió lineal és:

serie2_mall.lm <- lm(y~x, data= mall[1:51,])
summary(serie2_mall.lm)$adj.r.squared
## [1] 0.02816003

És a dir, un valor molt baix. Llavors aplicam la regressió segmentada seguint el mateix procediment que abans.

Vegem quants de punts de canvi tenim i identifiquem-los:

punts_canvi_serie2_mall <-selgmented(serie2_mall.lm, Kmax = 10, type="bic")
## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 7 .. 8 .. 9 .. 10 .. 
## 
## BIC to detect no. of breakpoints:
##        0        1        2     <NA>        3     <NA>        4        5 
## 489.3701 496.5603 503.9916 504.9916 505.6789 506.6789 500.9452 502.8412 
##        6        7        8 
## 496.1681 453.0991 440.6423 
## 
## No. of selected breakpoints: 8
serie2_mall.seg <- segmented(serie2_mall.lm, seg.Z=~x, psi = c(4,10,17,23,28,35,40,47))
summary(serie2_mall.seg)$psi
##        Initial      Est.    St.Err
## psi1.x       4  4.263718 0.7422268
## psi2.x      10 11.000003 0.5850117
## psi3.x      17 16.430363 0.4682119
## psi4.x      23 22.901107 0.4187605
## psi5.x      28 27.805715 0.4687236
## psi6.x      35 34.907009 0.4340051
## psi7.x      40 40.227949 0.3993340
## psi8.x      47 46.922405 0.3796441

Ens indica que tenim punts de canvi als mesos de 1-2016, 8-2016, 1-2017, 8-2017, 1-2018, 8-2018, 1-2019, 8-2019. És a dir segueixen el mateix patró que abans.

Ara el valor de \(R^2\) de la segmentada és

summary(serie2_mall.seg)$adj.r.squared
## [1] 0.8316885

Aquest és més proper a 1 i podem assegurar que és bo.

Anem a visualitzar la regressió segmentada damunt les nostres dades

p_serie2_mall <- ggplot(mall[1:51,], aes(x=x, y=y)) + geom_line()+#dibuixam les nostres dades en línia contínua + 
ggtitle('Regressió lineal i segmentada sobre les dades de Mallorca \ndurant la COVID') + xlab('índex del mes')+ ylab('Despeses mensuals en €')+ ylim(c(0,1300))

my.coef2_mall <- coef(serie2_mall.lm) #coeficients de la recta de regressió lineal

p_serie2_mall <- p_serie2_mall + geom_abline(intercept=my.coef2_mall[1], slope=my.coef2_mall[2], color="green") #afegim la recta de regressió lineal

my.model2_mall <- data.frame(x=1:51, y=fitted(serie2_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada

p_serie2_mall <- p_serie2_mall + geom_line(data=my.model2_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines2_mall <- serie2_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats

p_serie2_mall <- p_serie2_mall + geom_vline(xintercept=my.lines2_mall, linetype="dashed") #línies verticals en els punts de canvi

p_serie2_mall <- p_serie2_mall + geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

ggplotly(p_serie2_mall)

Per poder escriure la funció necessitam les pendents i els punts d’intersecció amb l’eix OY, que ens ho donen les següents funcions:

#PENDENTS
slope(serie2_mall.seg)
## $x
##           Est. St.Err. t value CI(95%).l CI(95%).u
## slope1 -38.707 21.3120 -1.8162   -82.065    4.6519
## slope2  40.282  9.0058  4.4729    21.959   58.6040
## slope3 -61.599 15.0700 -4.0877   -92.259  -30.9400
## slope4  62.878 11.3910  5.5197    39.702   86.0540
## slope5 -75.010 15.0700 -4.9776  -105.670  -44.3510
## slope6  49.277  9.0058  5.4717    30.954   67.5990
## slope7 -72.938 11.3910 -6.4028   -96.114  -49.7610
## slope8  66.757 11.3910  5.8603    43.581   89.9340
## slope9 -85.307 15.0700 -5.6609  -115.970  -54.6480
#PUNTS D'INTERSECCIÓ
intercept(serie2_mall.seg)
## $x
##                Est.
## intercept1   934.94
## intercept2   598.16
## intercept3  1718.80
## intercept4  -326.36
## intercept5  2831.40
## intercept6  -624.46
## intercept7  3641.70
## intercept8 -1978.00
## intercept9  5157.30

Aleshores tenim

\(\hat{y}= \left\{ \begin{array}{lcc} -38.707x + 934.94 & si & x \leq \psi_1 \\ \\ 40.282x + 598.16 & si & \psi_1 < x \leq \psi_2 \\ \\ -61.599x + 1718.8 & si & \psi_2 < x \leq \psi_3 \\ \\ 62.878x - 326.36 & si & \psi_3 < x \leq \psi_4 \\ \\ -75.010x + 2831.4 & si & \psi_4 < x \leq \psi_5 \\ \\ 49.277x - 624.46 & si & \psi_5 < x \leq \psi_6 \\ \\ -72.938x + 3641.7 & si & \psi_6 < x \leq \psi_7\\ \\ 66.757x - 1978 & si & \psi_7 < x \leq \psi_8\\ \\ -85.307x + 5157.3 & si & \psi_8 < x \\ \end{array} \right.\)

on \(\psi_1= 4.26\), \(\psi_2= 11\), \(\psi_3= 16.43\), \(\psi_4= 22.9\), \(\psi_5= 27.8\), \(\psi_6= 34.91\), \(\psi_7= 40.23\) i \(\psi_8= 46.92\).

L’error del model de regressió segmentada és (ens interessa RMSE):

accuracy(serie2_mall.seg)
##                        ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 4.458103e-15 38.33289 32.93155 -0.1657412 3.661674 0.3423069

És a dir 38.33289 euros.

Anem a fer la predicció, en aquest cas també tenim que els següents punts de canvi es donaran en gener, agost i gener.

  • La pendent de 8-2019 a 1-2020 i de 8-2020 a 1-2021 serà: (-61.599 - 75.010 -72.938)/3 = -69.849
  • La pendent de 1-2020 a 8-2020 serà : (40.282 + 62.878 + 49.277 + 66.757)/4 = 54.799

Seguint el mateix procediment que abans obtenim que \(n_1 = 4431.55\), \(n_2 = -2050.14\) i \(n_3 = 5304.09\).

És a dir la funció per a la predicció de 2020 és :

\(\hat{y}= \left\{ \begin{array}{lcc} -69.849x + 4431.55 & si & x \leq \psi_9 \\ \\ 54.799x - 2050.14 & si & \psi_9 < x \leq \psi_{10} \\ \\ -69.849x + 5304.09 & si & \psi_{10} < x \\ \end{array} \right.\)

on \(\psi_9 = 52\) i \(\psi_8 = 59\).

#Ho graficam

p2_serie_mall <- ggplot(mall, aes(x=x, y=y)) + geom_line()+  #dibuixam les nostres dades en línia contínua
  ggtitle('Predicció de Mallorca amb el model de regressió segmentada \ndurant la COVID') + 
  xlab('Índex del mes') + 
  ylab('Despeses mensuals en €') +
  ylim(c(0,1300)) #dibuixam les nostres dades en línia contínua

my.model2_mall <- data.frame(x=1:51, y=fitted(serie2_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada

p2_serie_mall <- p2_serie_mall + geom_line(data=my.model2_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines2_mall <- serie2_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats

p2_serie_mall <- p2_serie_mall + geom_vline(xintercept=my.lines2_mall, linetype="dashed") #línies verticals en els punts de canvi

p2_serie_mall <- p2_serie_mall + geom_vline(xintercept=0) + geom_hline(yintercept=0)+ theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

p2_serie_mall <- p2_serie_mall + geom_abline(intercept = 4431.55, slope=-69.849, colour="green") +                 
                                 geom_abline(intercept = -2050.14 , slope=54.799, colour="blue") +
                                 geom_abline(intercept = 5304.09 , slope=-69.849, colour="orange")

ggplotly(p2_serie_mall)

Podem observar gràficament que per a l’any 2020 la predicció no és gens bona ja que hi va haver un fet inesperat que és la COVID.

Calculem aquest error:

o2_mall<-c(serie_mall[52:63]) # dades reals per fer predicció del 2020

v3=c(52:59)
f3 <- sapply(v3, function(x) 54.799*x-2050.14) #predicció de gener 2020 a agost 2020  

v4=c(60:63)
f4 <- sapply(v4, function(x) -69.849*x + 5304.09) #predicció de setembre 2020 a desembre 2020

p2_mall <- c(f3,f4) #predicció de 2020

sqrt(sum((p2_mall-o2_mall)^2)/12)
## [1] 442.1241

En aquest cas l’error de la predicció amb aquest model és d’uns 425 euros.


- DECOMPOSE

Fem el mateix procediment que abans amb el mètode decompose i predim l’any 2020.

La descomposició, en aquest cas, és la següent:

decompose_s2_mall <- decompose(serie2_mall)
plot(decompose_s2_mall,xlab="Any")

On les components de tendència són:

t2_mall <- decompose_s2_mall$trend #tendències de la sèrie sense el tros a predir abans del COVID 
t2_mall
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2015                                                                        
## 2016       NA       NA       NA 899.8458 899.0708 896.9637 895.5496 888.9975
## 2017 891.8125 896.2783 901.1758 905.1896 907.7217 908.4975 909.7354 914.6108
## 2018 924.3500 924.6162 923.1429 921.8554 921.3479 920.0171 918.4050 914.4658
## 2019 914.9921 918.6238 923.0717 925.9079 928.4362 930.7558       NA       NA
##           Sep      Oct      Nov      Dec
## 2015                NA       NA       NA
## 2016 885.1633 886.3112 887.5750 888.9183
## 2017 920.9775 923.1212 922.9363 923.5758
## 2018 909.9675 910.4688 911.8138 912.6679
## 2019       NA       NA       NA       NA

I les estacionals:

s2_mall <- decompose_s2_mall$seasonal #estacionalitat
s2_mall <-  s2_mall[4:15] #estacionalitat de gener a desembre
s2_mall
##  [1] -174.403157 -122.557740  -54.608435  -60.967983    3.885038   51.915663
##  [7]  179.218371  198.616982   82.055593   71.194621 -116.073296  -58.275657

Fem la predicció:

a2_mall <- c(s2_mall[7:12],s2_mall) #estacionalitat de juliol 2019 a desembre 2020 (es per poder fer la predicció)

pred2_decompose_mall <- sapply(a2_mall, function(x) 930.7558 + x) #predicció de juliol de 2019 a desembre 2020

p2_dec_mall<-c(serie2_mall[1:45], pred2_decompose_mall)

A2_mall<- data.frame("x" = mall[1:63,]$x, "y" = mall[1:63,]$y, "p"= p2_dec_mall)

grafica2_mall_dec <- ggplot(A2_mall)+
  geom_line(aes(x,p), color="red")+
  geom_line(aes(x,y))+
  labs(title="Predicció de Mallorca durant la COVID amb el model de descomposició", x="Índex del mes", y="Despesese mensuals en €")+
  geom_vline(xintercept = 0)+ geom_hline(yintercept = 0)+ theme(panel.background = element_blank())

grafica2_mall_dec

Amb aquest mètode també s’observa que la caiguda del COVID no s’ha detectat, ja que de les dades on estam aprenent per fer la predicció no hi tenim valors pareguts als de la COVID.

L’error que es comet és:

sqrt(sum((c(serie_mall[52:63]- pred2_decompose_mall[7:18]))^2)/12)
## [1] 390.5974

Que també és un valor gros, tot i que no tant com amb el model anterior.


Vegem, igual que abans, la predicció amb la funció predict():

prediccio2_mall <- predict(serie2_mall,h=12)
plot(prediccio2_mall, xlab="Any", ylab ="Despeses menusals en €")

Així pareix que ha de ser una bona predicció, però recordem que aquest cicle és el de la COVID, vegem com queda la predicció i la sèrie original:

df_prediccio2_mall <- data.frame(prediccio2_mall)
p2_ets_mall <- data.frame("x"= 1:63, "PointForecast"=c(serie2_mall,df_prediccio2_mall$Point.Forecast), "Lo80" = c(rep(NA,51),df_prediccio2_mall$Lo.80), "Hi80" = c(rep(NA,51),df_prediccio2_mall$Hi.80), "Lo95" = c(rep(NA,51),df_prediccio2_mall$Lo.95),"Hi95" = c(rep(NA,51),df_prediccio2_mall$Hi.95))

grafica_pred2_ets <- ggplot((mall[1:63,]))+
  geom_ribbon(data = p2_ets_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data =p2_ets_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = p2_ets_mall, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció del model ETS(A,N,A) a Mallorca durant la COVID", y="Despeses mensuals en €", x="Índex del mes") + 
  geom_vline(xintercept = 0) + geom_hline(yintercept = 0)+ theme(panel.background = element_blank())

grafica_pred2_ets

Efectivament és una mala predicció, ja que no ha detectat la baixada de la COVID.

Si calculam l’error comés, aquest és de quasi 400 euros.

accuracy(prediccio2_mall,serie3_mall)
##                       ME      RMSE       MAE       MPE     MAPE      MASE
## Training set    1.768898  27.38848  22.10859 0.1014835 2.474142 0.6564259
## Test set     -248.774469 392.70874 248.77447      -Inf      Inf 7.3863591
##                     ACF1 Theil's U
## Training set 0.004436559        NA
## Test set     0.439341466         0


- SARIMA

Quin model proposam nosaltres? Vegem l’ACF i el PACF.

par(mfrow=c(1,2))
acf(serie2_mall)
pacf(serie2_mall)

Per a la part regular, de l’ACF obtenim un MA(2) i del PACF un AR(1). Notem que no fa falta fer cap diferenciació a la part regular ja que no considram que hi hagi tendència. No obstant, sí que observam que hi ha una certa estacionalitat per la forma que tenen les barres en els correlogrames.

Llavors fem una diferenciació d’ordre 12:

serie2_mall_diff <- diff(serie2_mall,12)
plot(serie2_mall_diff, main="Sèrie sense estacionalitat", xlab="Any", ylab="Sèrie diferenciada")

Ara ja no s’observa l’estacionalitat, llavors hem de fer una diferenciació D=1. Vegem, ara, quins són els correlogrames

par(mfrow=c(1,2))
acf(serie2_mall_diff,lag=36)
pacf(serie2_mall_diff,lag=36)

Amb aquests obtenim un MA(1) i un AR(1) per a la part estacional. És a dir, per predir un més hem de tenir en compte el valor i els errors comesos d’aquest mateix mes un cicle anterior.

Llavors el model que nosaltres proposam es un ARIMA(1,0,2)(1,1,1)[12]

R proposa el següent model:

auto.arima(serie2_mall)
## Series: serie2_mall 
## ARIMA(0,0,0)(1,1,0)[12] with drift 
## 
## Coefficients:
##          sar1   drift
##       -0.7165  0.9022
## s.e.   0.1136  0.2674
## 
## sigma^2 = 927:  log likelihood = -191.85
## AIC=389.7   AICc=390.39   BIC=394.69

Fixem-nos que és el mateix model que proposa R per a la serie1, és a dir, la que va d’octubre de 2015 a desembre de 2018.

Els dos models que tenim són:

model3_mall <- arima(serie2_mall, order=c(1,0,2), seasonal = list(order=c(1,1,1), period=12))

model4_mall <- arima(serie2_mall, order=c(0,0,0), seasonal = list( order=c(1,1,0), period=12))

I els seus errors

accuracy(model3_mall)
##                    ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 5.895861 22.21249 15.13338 0.5278811 1.683231 0.2080521
##                     ACF1
## Training set -0.02958489
accuracy(model4_mall)
##                    ME     RMSE      MAE       MPE     MAPE      MASE      ACF1
## Training set 10.85733 30.49463 21.65555 0.9581157 2.367599 0.2977183 0.2066374

Les prediccions són:

fc_m3_mall <- forecast(model3_mall, h=12)
fc_m4_mall <- forecast(model4_mall,h=12)
#PRIMER MODEL
pre2_mall <- data.frame("Point Forecast" = serie2_mall, "Lo 80" =  rep(NA,51), "Hi 80"= rep(NA,51), "Lo 95" = rep(NA,51), "Hi 95" = rep(NA,51))

pred_m3_mall <-data.frame(fc_m3_mall)

grafica_m3_mall <- data.frame("x" = 1:63, "PointForecast" = c(pre2_mall$Point.Forecast,pred_m3_mall$Point.Forecast), "Lo80" = c(pre2_mall$Lo.80, pred_m3_mall$Lo.80), "Hi80"= c(pre2_mall$Hi.80, pred_m3_mall$Hi.80), "Lo95" = c(pre2_mall$Lo.95, pred_m3_mall$Lo.95), "Hi95" = c(pre2_mall$Hi.95, pred_m3_mall$Hi.95))


grafica3_mall <- ggplot(mall[1:63,])+
  geom_ribbon(data = grafica_m3_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m3_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m3_mall, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red") +
  labs(title="Prediccó  de Mallorca amb el model ARIMA (1,0,2)(1,1,1)[12] durant la COVID", x="Índex del mes", y="Despeses mensuals en €") +
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica3_mall

#SEGON MODEL
pred_m4_mall <-data.frame(fc_m4_mall)

grafica_m4_mall <- data.frame("x" = 1:63, "PointForecast" = c(pre2_mall$Point.Forecast,pred_m4_mall$Point.Forecast), "Lo80" = c(pre2_mall$Lo.80, pred_m4_mall$Lo.80), "Hi80"= c(pre2_mall$Hi.80, pred_m4_mall$Hi.80), "Lo95" = c(pre2_mall$Lo.95, pred_m4_mall$Lo.95), "Hi95" = c(pre2_mall$Hi.95, pred_m4_mall$Hi.95))

grafica4_mall <- ggplot(mall[1:63,])+
  geom_ribbon(data = grafica_m4_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m4_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m4_mall, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red") +
  labs(title="Predicció de Mallorca amb el model ARIMA (0,0,0)(1,1,0)[12] durant la COVID", x="Índex del mes", y="Despeses mensuals en €") +
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica4_mall

Amb un error de 384.5 i 383.6:

accuracy(fc_m3_mall,serie_mall[52:63], h=12)
##                       ME      RMSE       MAE       MPE     MAPE      MASE
## Training set    5.895861  22.21249  15.13338 0.5278811 1.683231 0.2080521
## Test set     -243.110836 384.52595 243.11084      -Inf      Inf 3.3422626
##                     ACF1
## Training set -0.02958489
## Test set              NA
accuracy(fc_m4_mall,serie_mall[52:63], h=12)
##                      ME      RMSE       MAE       MPE     MAPE      MASE
## Training set   10.85733  30.49463  21.65555 0.9581157 2.367599 0.2977183
## Test set     -237.63235 383.58017 238.11170      -Inf      Inf 3.2735350
##                   ACF1
## Training set 0.2066374
## Test set            NA

En els dos models veiem que la predicció és dolenta ja que segueixen el mateix patró que els cicles anteriors i no preveuen la COVID. També podem observar que l’error que es comet és bastant alt.

Igual que abans, esdudiem amb més detall aquests errors:

checkresiduals(fc_m3_mall)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,2)(1,1,1)[12]
## Q* = 4.7827, df = 5, p-value = 0.443
## 
## Model df: 5.   Total lags used: 10
checkresiduals(fc_m4_mall)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(0,0,0)(1,1,0)[12]
## Q* = 13.035, df = 9, p-value = 0.161
## 
## Model df: 1.   Total lags used: 10
par(mfrow=c(1,2))
qqPlot(fc_m3_mall$residuals, id=FALSE, mean=mean(fc_m3_mall$residuals), sd=sd(fc_m3_mall$residuals))
qqPlot(fc_m4_mall$residuals, id=FALSE, mean=mean(fc_m4_mall$residuals), sd=sd(fc_m4_mall$residuals))

shapiro.test(fc_m3_mall$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m3_mall$residuals
## W = 0.93903, p-value = 0.01123
shapiro.test(fc_m4_mall$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m4_mall$residuals
## W = 0.93172, p-value = 0.005797
lillie.test(fc_m3_mall$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m3_mall$residuals
## D = 0.14786, p-value = 0.00709
lillie.test(fc_m4_mall$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m4_mall$residuals
## D = 0.20836, p-value = 7.9e-06

Obtenim que els tests de normalitat ens donen uns p-valors molt petits. No obstant, en els dos models els residus no segueixen cap estructura i les barres dels ACF es mantenen dins les bandes de confiança. També observam que, als qqPlot, tot i haver-hi valors fora de la regió blava i de no estar sobre la diagonal, la majoria de valors es troben al voltant de 0.


Vegem un resum dels errors que cometen cada un dels models en aquest cas:

reg. segmentada descomposició clàssica ETS(A,N,A) ARIMA (1,0,2)(1,1,1)[12] ARIMA (0,0,0)(1,1,0)[12]
error model 38.33 27.39 22.21 30.49
error predicció 424.78 390.6 392.71 384.53 383.58


2.3. DESPRÉS DE LA COVID


- REGRESSIÓ SEGMENTADA

El valor de \(R^2\) per a la regressió lineal és

serie3_mall.lm <- lm(y~x, data= mall[1:63,])
summary(serie3_mall.lm)$adj.r.squared
## [1] 0.02681106

Vegem els punts de canvi:

punts_canvi_serie3_mall <-selgmented(serie3_mall.lm, Kmax = 10, type="bic")
## No. of breakpoints: 2 .. 3 .. 4 .. 5 .. 6 .. 
## (search truncated at 6 breakpoints due to increasing values of BIC) 
## 
## BIC to detect no. of breakpoints:
##        0        1        2        3        4        5        6 
## 670.4792 673.4233 681.6172 688.4189 695.2176 701.0818 706.7486 
## 
## No. of selected breakpoints:  0

Fixem-nos que si demanam a R que ens indiqui quants de punts de canvi estima amb la sèrie completa, és a dir, amb la COVID, no detecta cap punt de canvi. Aquest fet s’eslica amb detall al treball.

serie3_mall.seg <- segmented(serie3_mall.lm, seg.Z=~x, psi = c(5,10,17,23,28,35,40,47,55,58))
summary(serie3_mall.seg)$psi
##         Initial      Est.    St.Err
## psi1.x        5  5.783722 1.7320089
## psi2.x       10 10.400505 1.1467103
## psi3.x       17 16.853561 1.0053175
## psi4.x       23 22.744051 0.8995765
## psi5.x       28 27.944118 1.0394421
## psi6.x       35 34.740660 0.9251495
## psi7.x       40 40.250732 0.7944349
## psi8.x       47 47.223848 0.5600066
## psi9.x       55 56.000020 0.4535518
## psi10.x      58 58.017477 0.3070058

Aquests es donen als mesos següents: 3-2026, 7-2016, 2-2017, 8-2017, 1-2018, 8-2018, 1-2019, 8-2019, 5-2020, 7-2020.

Ara, el valor de l’\(R^2\) de la segmentada és

summary(serie3_mall.seg)$adj.r.squared
## [1] 0.7222148

Anem a visualitzar la regressió segmentada damunt les nostres dades

p_serie3_mall <- ggplot(mall[1:63,], aes(x=x, y=y)) + geom_line()+ #dibuixam les nostres dades en línia contínua
  ggtitle('Regressió lineal i segmentada sobre les dades \nde Mallorca després de la COVID') + xlab('índex del mes')+ ylab('Despeses mensuals en €')+
  ylim(c(0,1300))
  #dibuixam les nostres dades en línia contínua

my.coef3_mall <- coef(serie3_mall.lm) #coeficients de la recta de regressió lineal

p_serie3_mall <- p_serie3_mall + geom_abline(intercept=my.coef3_mall[1], slope=my.coef3_mall[2], color="green") #afegim la recta de regressió lineal

my.model3_mall <- data.frame(x=1:63, y=fitted(serie3_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada

p_serie3_mall <- p_serie3_mall + geom_line(data=my.model3_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines3_mall <- serie3_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats

p_serie3_mall <- p_serie3_mall + geom_vline(xintercept=my.lines3_mall, linetype="dashed") #línies verticals en els punts de canvi

p_serie3_mall <- p_serie3_mall + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0) + theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic 

ggplotly(p_serie3_mall)

Per escriure la funció a trossos tenim:

#PUNTS D'INTERSECCIÓ
intercept(serie3_mall.seg)
## $x
##                  Est.
## intercept1     897.69
## intercept2     449.26
## intercept3    1601.80
## intercept4    -444.27
## intercept5    2752.50
## intercept6    -722.00
## intercept7    3636.00
## intercept8   -2243.00
## intercept9    6553.70
## intercept10 -23235.00
## intercept11   5426.40
#PENDENTS
slope(serie3_mall.seg)
## $x
##             Est. St.Err.  t value CI(95%).l CI(95%).u
## slope1   -19.420  32.900 -0.59027   -85.863   47.0230
## slope2    58.115  32.900  1.76640    -8.328  124.5600
## slope3   -52.703  24.870 -2.11910  -102.930   -2.4764
## slope4    68.701  24.870  2.76240    18.475  118.9300
## slope5   -71.853  32.900 -2.18400  -138.300   -5.4096
## slope6    52.484  19.662  2.66940    12.777   92.1920
## slope7   -72.961  24.870 -2.93370  -123.190  -22.7340
## slope8    73.101  19.662  3.71800    33.394  112.8100
## slope9  -113.180  13.431 -8.42630  -140.300  -86.0520
## slope10  418.770 147.130  2.84620   121.630  715.9100
## slope11  -75.248  32.900 -2.28720  -141.690   -8.8050

L’error de la regressió és:

accuracy(serie3_mall.seg)
##                        ME     RMSE     MAE  MPE MAPE      MASE
## Training set 3.609932e-15 83.93007 53.6569 -Inf  Inf 0.4360175

Anem a fer la predicció per a l’any 2021. Recordem que els punts de canvi es donen el 3-2026, 7-2016, 2-2017, 8-2017, 1-2018, 8-2018, 1-2019, 8-2019, 5-2020, 7-2020. Degut a la pertorbació del COVID, sí que hi segueix havent un punt de canvi en estiu i l’altre en hivern successivament però ara no és donen al mateix mes. Per això, per predir els següents punts de canvi agafarem el mes més freqüent. Aleshores els següents punts de canvi seran en 1-2021, 8-2021 i 1-2022.

  • La pendent de 7-2020 a 1-2021 i de 8-2021 a 1-2022 serà : -77.674
  • La pendent de 1-2021 a 8-2021 serà : 134.234

Els nous punts d’intersecció es calculen de forma anàloga que als casos anteriors. Per a les tres noves rectes obtenim que són \(n_1=5563.7164\), \(n_2=-7998.396\) i \(n_3=7047.07\).

Aleshores, la predicció de 2021 serà:

\(\hat{y}= \left\{ \begin{array}{lcc} -77.674x + 5563.716 & si & x \leq \psi_{11} \\ \\ 134.234x - 7998.396 & si & \psi_{11} < x \leq \psi_{12} \\ \\ -77.674x + 7047.07 & si & \psi_{12} < x \\ \end{array} \right.\)

on \(\psi_7 = 64\) i \(\psi_8 = 71\).

#visualitzem-la 

p3_serie_mall <- ggplot(mall, aes(x=x, y=y)) + geom_line()+ #dibuixam les nostres dades en línia contínua
  ggtitle('Predicció de Mallorca amb el model de \nregressió segmentada després de la COVID') + xlab('índex del mes')+ ylab('Despeses mensuals en €')+
  ylim(c(0,1600))

my.model3_mall <- data.frame(x=1:63, y=fitted(serie3_mall.seg)) #model nou amb els valors ajustats de la regressió segmentada

p3_serie_mall <- p3_serie_mall + geom_line(data=my.model3_mall, aes(x=x,y=y), colour="red") #afegim la recta de regressió segmentada

my.lines3_mall <- serie3_mall.seg$psi[,2]# guardam on estan els punts de canvi estimats

p3_serie_mall <- p3_serie_mall + geom_vline(xintercept=my.lines3_mall, linetype="dashed") #línies verticals en els punts de canvi

p3_serie_mall <- p3_serie_mall + geom_vline(xintercept = 0)+ geom_hline(yintercept = 0)+theme(panel.background = element_blank()) #eliminam el fons i la quadrícula del gràfic

p3_serie_mall <- p3_serie_mall + geom_abline(intercept = 5563.7164, slope=-77.674, colour="green") +                 
                                 geom_abline(intercept = -7998.396 , slope= 134.234, colour="blue") +
                                 geom_abline(intercept = 7047.07 , slope=-77.674, colour="orange")

ggplotly(p3_serie_mall)

A simple vista observam que la predicció no és massa bona, ja que la predicció que es prediu pel mes d’estiu de l’any 2021 és d’uns 1500€ aproximadament, mentre que el valor real és d’uns 1100€.

Vegem quin és aquest error que es comet:

o3_mall<-c(serie_mall[64:75]) # dades reals per fer predicció del 2021

v5=c(64:71)
f5 <- sapply(v5, function(x) 134.234*x-7998.396) #predicció de gener 2021 a agost 2021 

v6=c(72:75)
f6 <- sapply(v6, function(x) -77.674*x + 7047.07) #predicció de setembre 2021 a desembre 2021


p3_mall <- c(f5,f6) #predicció de 2021

sqrt(sum((p3_mall-o3_mall)^2)/12)
## [1] 269.9865


- DECOMPOSE

Visualitzem la sèrie descomposada:

decompose_s3_mall <- decompose(serie3_mall)
plot(decompose_s3_mall, xlab="Any")

Els components de tendència són:

t3_mall <- decompose_s3_mall$trend #tendències de la sèrie sense el tros a predir abans de la COVID
t3_mall
##           Jan      Feb      Mar      Apr      May      Jun      Jul      Aug
## 2015                                                                        
## 2016       NA       NA       NA 899.8458 899.0708 896.9637 895.5496 888.9975
## 2017 891.8125 896.2783 901.1758 905.1896 907.7217 908.4975 909.7354 914.6108
## 2018 924.3500 924.6162 923.1429 921.8554 921.3479 920.0171 918.4050 914.4658
## 2019 914.9921 918.6238 923.0717 925.9079 928.4362 930.7558 931.6700 934.0650
## 2020 750.7167 737.9196 727.9308 713.3075 700.2629 692.3146       NA       NA
##           Sep      Oct      Nov      Dec
## 2015                NA       NA       NA
## 2016 885.1633 886.3112 887.5750 888.9183
## 2017 920.9775 923.1212 922.9363 923.5758
## 2018 909.9675 910.4688 911.8138 912.6679
## 2019 931.8875 891.1542 816.3763 767.6412
## 2020       NA       NA       NA       NA

I els d’estacionalitat:

s3_mall <- decompose_s3_mall$seasonal #estacionalitat
s3_mall <-  s3_mall[4:15] #estacionalitat de gener a desembre
s3_mall
##  [1] -127.84021  -68.39937  -27.68521 -186.60614 -132.11481   55.96136
##  [7]  191.44261  213.85281   92.39615   86.61625  -78.02771  -19.59573

Anem a fer la predicció:

a3_mall <- c(s3_mall[7:12],s3_mall) #estacionalitat de juliol 2020 a desembre 2021 (es per poder fer la predicció)

pred3_decompose_mall <- sapply(a3_mall, function(x) 692.3146 + x) #predicció de juliol de 2019 a desembre 2020

p3_dec_mall<-c(serie3_mall[1:57], pred3_decompose_mall)

A3_mall<- data.frame("x" = mall[1:75,]$x, "y" = mall[1:75,]$y, "p"= p3_dec_mall)

grafica3_mall_dec <- ggplot(A3_mall)+
  geom_line(aes(x,p), color="red")+
  geom_line(aes(x,y))+
  labs(title="Predicció de Mallorca després de la COVID amb el model de descomposició", x="Índex del mes", y="Despeses mensualns en €")+geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica3_mall_dec

Observam que la predicció, en vermell, segueix més o menys el mateix patró que la sèrie original. No obstant això, es troba a un nivell més baix degut a les dades anteriors que té per aprendre, on hi ha l’efecte de la COVID.

Calculem l’error de la predicció:

sqrt(sum((c(serie_mall[64:75]- pred3_decompose_mall[7:18]))^2)/12)
## [1] 280.1672


Amb la funció predict() la predicció seria la següent:

prediccio3_mall <- predict(serie3_mall,h=12)
plot(prediccio3_mall, xlab="Any", ylab="Despeses mensuals en €")

Observam que, a diferència de les series anteriors amb aquesta funció, sense la COVID, la predicció és molt dolenta. Vegem-la juntament amb les nostres dades:

df_prediccio3_mall <- data.frame(prediccio3_mall)
p3_ets_mall <- data.frame("x"= 1:75, "PointForecast"=c(serie3_mall,df_prediccio3_mall$Point.Forecast), "Lo80" = c(rep(NA,63),df_prediccio3_mall$Lo.80), "Hi80" = c(rep(NA,63),df_prediccio3_mall$Hi.80), "Lo95" = c(rep(NA,63),df_prediccio3_mall$Lo.95),"Hi95" = c(rep(NA,63),df_prediccio3_mall$Hi.95))

grafica_pred3_ets <- ggplot((mall[1:75,]))+
  geom_ribbon(data = p3_ets_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data =p3_ets_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = p3_ets_mall, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció del model ETS(A,N,N) a Mallorca després de la COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica_pred3_ets

Aquest fet és degut a que el model ETS dóna més pes a les observacions anteriors properes i no a les llunyanes, i no està considerant estacionalitat degut a la detecció de la COVID.

accuracy(prediccio3_mall, serie_mall)
##                      ME     RMSE       MAE      MPE     MAPE     MASE      ACF1
## Training set  -3.019093 162.3973  95.90263     -Inf      Inf 1.119973 0.1219558
## Test set     251.816426 282.9966 251.81643 24.85994 24.85994 2.940770 0.6915062
##              Theil's U
## Training set        NA
## Test set      3.250169


- SARIMA

Quin model proposam nosaltres? Vegem l’ACF i el PACF:

par(mfrow=c(1,2))
acf(serie3_mall, lag.max=36)
pacf(serie3_mall, lag.max=36)

Per a la part regular obtenim un MA(1), un AR(2). Observam una clara estacionalitat el l’ACF però com les barres no surten no la consideram significativa. Llavors d=1.

Diferenciem la part estacional i calculem el seu ACF i PACF:

serie3_mall_diff <- diff(serie3_mall,12)
plot.ts(serie3_mall_diff, main="Sèrie sense estacionalitat", ylab="Sèrie diferenciada", xlab="Any")

par(mfrow=c(1,2))
acf(serie3_mall_diff)
pacf(serie3_mall_diff)

Fent una diferenciació d’ordre 12, per llevar estacionalitat, i tornant a calcular els ACF i PACF obtenim: que l’ordre del model autoregressiu és P=2 i el de mitjana mòbil Q=1.

Aleshores la nostra proposta és un model ARIMA(2,1,1)(2,0,1)[12]

R proposa el següent model:

auto.arima(serie3_mall)
## Series: serie3_mall 
## ARIMA(1,0,1)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          ar1     ma1    sar1      mean
##       0.3171  0.7225  0.6045  831.5038
## s.e.  0.1369  0.0994  0.1267   75.8692
## 
## sigma^2 = 14934:  log likelihood = -393.43
## AIC=796.86   AICc=797.91   BIC=807.58

Llavors tenim els següents models:

model5_mall <- arima(serie3_mall, order=c(2,1,1), seasonal = list( order=c(2,0,1), period=12))

model6_mall <- arima(serie3_mall, order=c(1,0,1), seasonal = list( order=c(1,0,0), period=12))

Amb un error de 112.18 i 118.26, respectivament:

accuracy(model5_mall)
##                     ME     RMSE      MAE MPE MAPE      MASE        ACF1
## Training set -3.612028 127.6573 76.78002 NaN  Inf 0.7880336 -0.03300271
accuracy(model6_mall)
##                     ME     RMSE      MAE  MPE MAPE      MASE         ACF1
## Training set -4.030176 118.2628 71.00525 -Inf  Inf 0.7287641 0.0005558812

Les prediccions d’aquests models són:

fc_m5_mall <- forecast(model5_mall, h=12)

fc_m6_mall <- forecast(model6_mall,h=12)
#grafiquem-les

#primer model

pre3_mall <- data.frame("Point Forecast" = serie3_mall, "Lo 80" =  rep(NA,63), "Hi 80"= rep(NA,63), "Lo 95" = rep(NA,63), "Hi 95" = rep(NA,63))

pred_m5_mall <-data.frame(fc_m5_mall)


grafica_m5_mall <- data.frame("x" = 1:75, "PointForecast" = c(pre3_mall$Point.Forecast,pred_m5_mall$Point.Forecast), "Lo80" = c(pre3_mall$Lo.80, pred_m5_mall$Lo.80), "Hi80"= c(pre3_mall$Hi.80, pred_m5_mall$Hi.80), "Lo95" = c(pre3_mall$Lo.95, pred_m5_mall$Lo.95), "Hi95" = c(pre3_mall$Hi.95, pred_m5_mall$Hi.95))


grafica5_mall <- ggplot(mall[1:75,])+
  geom_ribbon(data = grafica_m5_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m5_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m5_mall, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció de Mallorca amb el model ARIMA (2,1,1)(2,0,1)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €")+ 
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica5_mall

#segon model 

pred_m6_mall <-data.frame(fc_m6_mall)

grafica_m6_mall <- data.frame("x" = 1:75, "PointForecast" = c(pre3_mall$Point.Forecast,pred_m6_mall$Point.Forecast), "Lo80" = c(pre3_mall$Lo.80, pred_m6_mall$Lo.80), "Hi80"= c(pre3_mall$Hi.80, pred_m6_mall$Hi.80), "Lo95" = c(pre3_mall$Lo.95, pred_m6_mall$Lo.95), "Hi95" = c(pre3_mall$Hi.95, pred_m6_mall$Hi.95))

grafica6_mall <- ggplot(mall[1:75,])+
  geom_ribbon(data = grafica_m6_mall, aes(x, ymin = Lo95, ymax = Hi95), fill = "blue", alpha = 0.25) +
  geom_ribbon(data = grafica_m6_mall, aes(x, ymin = Lo80, ymax = Hi80), fill = "blue", alpha = 0.25) +
  geom_line(data = grafica_m6_mall, aes(x, PointForecast), colour = "blue") +
  geom_line(aes(x,y), color="red")+
  labs(title="Predicció de Mallorca amb el model ARIMA (1,0,1)(1,0,0)[12] \ndesprés de la COVID", x="Índex del mes", y="Despeses mensuals en €") + 
  geom_vline(xintercept=0) + geom_hline(yintercept=0) + theme(panel.background = element_blank())

grafica6_mall

Tant a d’un model com a l’altre observam que la predicció no és molt bona. De la mateixa manera que als models anteriors, aquest fet és degut a l’alteració de les dades per la COVID.

Vegem quin és l’error de la predicció per a cada model

accuracy(fc_m5_mall,serie_mall[64:75], h=12)
##                      ME     RMSE       MAE      MPE     MAPE      MASE
## Training set  -3.612028 127.6573  76.78002      NaN      Inf 0.7880336
## Test set     422.131202 475.4906 422.13120 43.52191 43.52191 4.3325538
##                     ACF1
## Training set -0.03300271
## Test set              NA
accuracy(fc_m6_mall,serie_mall[64:75], h=12)
##                      ME     RMSE       MAE      MPE     MAPE      MASE
## Training set  -4.030176 118.2628  71.00525     -Inf      Inf 0.7287641
## Test set     236.006237 308.5205 240.53442 24.01452 24.63381 2.4687309
##                      ACF1
## Training set 0.0005558812
## Test set               NA

Obtenim que l’error és de 475.5 i 308.5.

Estudiem-los:

checkresiduals(fc_m5_mall)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,1,1)(2,0,1)[12]
## Q* = 9.1059, df = 7, p-value = 0.2451
## 
## Model df: 6.   Total lags used: 13
checkresiduals(fc_m6_mall)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(1,0,1)(1,0,0)[12] with non-zero mean
## Q* = 4.3222, df = 10, p-value = 0.9316
## 
## Model df: 3.   Total lags used: 13
par(mfrow=c(1,2))
qqPlot(fc_m5_mall$residuals, id=FALSE, mean=mean(fc_m5_mall$residuals), sd=sd(fc_m5_mall$residuals))
qqPlot(fc_m6_mall$residuals, id=FALSE, mean=mean(fc_m6_mall$residuals), sd=sd(fc_m6_mall$residuals))

shapiro.test(fc_m5_mall$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m5_mall$residuals
## W = 0.79747, p-value = 6.838e-08
shapiro.test(fc_m6_mall$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  fc_m6_mall$residuals
## W = 0.69775, p-value = 4.259e-10
lillie.test(fc_m5_mall$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m5_mall$residuals
## D = 0.18624, p-value = 1.068e-05
lillie.test(fc_m6_mall$residuals)
## 
##  Lilliefors (Kolmogorov-Smirnov) normality test
## 
## data:  fc_m6_mall$residuals
## D = 0.16427, p-value = 0.0002185

De la mateixa manera que als casos anteriors, els errors no segueixen una normal.


Vegem un resum dels errors de tots els models per aquest cas:

reg. segmentada descomposició clàssica ETS(A,N,N) ARIMA (2,1,1)(2,0,1)[12] ARIMA (1,0,1)(1,0,0)[12]
error model 83.93 162.4 112.18 118.26
error predicció 269.987 280.17 283 451.8 308.52